%clear all;
warning('off', 'MATLAB:colon:nonIntegerIndex')

%% segement all the files
FileList = {
   
  %  'mohEEG111.mat'
    'Marion12_2504.mat'
  %  'Tim12_2804.mat'
    %  'Romain12_2904.mat'
   % 'david12-2904.mat'
    };

%FileList = { 'session2.acq.mat' };
%Fs = 2000; default
Fs = 2000;

window_learn = 3000; % 3000 = 1.5 seconds
window_test = 1000;
shift = 0; % 1200 or 0 for marion
shift_mv_test = 0;
rest_window = 6000;
threshold = 8.4;% 8.4

X = []; % contains the trials from all classes
Y = []; % contains the classification from all classes
X_rest = [];
Y_rest = [];

is_pulsar_single_trigger_channel_from_eeglab = false;

for curFileIx = 1:length(FileList)
    

        load(['D:\temp\subjects\' FileList{curFileIx}]);
        disp(['Subject file:' FileList{curFileIx}]);

        datar = data(:,1:size(data,2)-2); % remove 2 trigger channel at the end

        %trial_channel = data(:,8);
        trial_channel = data(:, size(data,2)-1);

        %class_separation_channel = data(:,9);
        class_separation_channel = data(:, size(data,2));

        classN = length(count_sequence(class_separation_channel,5));  % "5" tags each class


        classRest = classN; % the rest class is always the last

        trialsN = length(count_sequence(trial_channel,5));  % "5" tags each trial
            
        disp(['Classes/Movements detected: ' int2str(classN)]);
        disp(['Total trails detected: ' int2str(trialsN)]);
        disp(['Frequency: ' int2str(Fs)]);

        for m = 1:classN
            
            %Pulsar BIOPAC 
            [best_index,best_trigger] = generate_trials_postitions_using_discrete_channel(m, trial_channel,class_separation_channel);
            
            % cut signal
            %index = find(diff(trigger)==1); % because we do it earlier now!!!!
            trials_per_class = zeros(size(datar,2),window_learn,length(best_index));
            
            trials_per_class_rest = zeros(size(datar,2),rest_window,length(best_index));
            disp(['Trials: ' int2str(length(best_index)) ' ' FileList{curFileIx}])

            %get the data for these triggers
            for i=1:length(best_index)
                
                %calculate a range around the trigger
                if (m == classRest)
                    buffer = datar(best_index(i):best_index(i)+rest_window-1,:);
                    trials_per_class_rest(:,:,i) = buffer';
                else
                    buffer = datar(best_index(i)+shift:best_index(i)+shift+window_learn-1,:);
                    trials_per_class(:,:,i) = buffer';
                end;
            end

            if (m == classRest)
               X_rest = cat(3,X_rest,trials_per_class_rest);
               Y_rest = cat(1,Y_rest,ones(length(best_index),1) * m);
            else
               X = cat(3,X,trials_per_class); % put all trials together
               Y = cat(1,Y,ones(length(best_index),1) * m);
            end;
            

        end;
   
    disp('--------------------------------');% next file
end

%% Train on X in CallClass

% the standard MDM

method_mean = 'riemann';
method_dist = 'riemann';

COVtrain = covariances(X);
Ytrain = Y;
labels = unique(Ytrain);
Nclass = length(labels);
CallClass = cell(Nclass,1);

% estimation of center - actual classification
for i=1:Nclass
    CallClass{i} = mean_covariances(COVtrain(:,:,Ytrain==labels(i)),method_mean);
end

%% Get test data

% we train a potato to detect movements
overlap = window_test-window_test/8;

%COV_all_signal=eeg2cov(datar',window,overlap);
%COV_movement = covariances(X);
COV_rest = covariances(X_rest);

%[C_mv, th] = potato_estimation_iterativ(COV_movement); % train on movement

[C_rest, th_rest] = potato_estimation_iterativ(COV_rest); % train on rest

[COVtest CovTestIndex] = eeg2cov_index(datar',window_test,overlap);

%[~,d_mv] = potato_detection(COVtest,C_mv,th);

[~,d_rest] = potato_detection(COVtest,C_rest,th_rest);


trigger = d_rest>threshold;

current_index = find(diff(trigger)==1);

trials = length(current_index);

disp(['Test trials detected: ' int2str(trials) ]);
%% detect the good threshold for movement detection
% n_trials = (classN-1) * 20;
% 
% last = 10000;
% best = -1;
% 
% for i = 1.40:0.1:8
%    
%    %skip if the fastest movement is less than ... 400 ms = 800 samples
%    minSamples = 600; % 300ms
%    minn = window / minSamples; % 1 in d equals 1 window, so we convert minn to windows too
%    
%    %if the smallest detected d is less than a movement, then skip all
%    if min(count_sequence(d_rest>i,1)) < minn 
%        continue;
%    end;
%    
%    if (abs((length(count_sequence(d_rest>i,1)) - n_trials)) < last)
%        last = abs(length(count_sequence(d_rest>i,1)) - n_trials);
%        best = i
%    end;
% 
% end;
% 
% th = best;
% %th = th_rest;
% 
% NTesttrial = count_sequence(d_rest>th,1);%sum(d_rest>th)
%  
% COV = zeros(10,10,NTesttrial);
% 
% % now we use the potato to extract movements from our signal 
pos = 0;
COV = [];
%range_to_be_processed = d_rest>th;
 
for i = current_index'
    
    %if (range_to_be_processed(i)==1)
       
        buffer = datar(CovTestIndex(i) + shift_mv_test : CovTestIndex(i) + window_test + shift_mv_test,:);
        cc = covariances(buffer');
        pos = pos + 1;
        COV(:,:,pos) = cc;
    %end;
 
end;

%% Now we have movements detected and trained classifier in CallClass

COVtest = COV;

d_test = zeros(length(current_index),Nclass);
for j=1:length(current_index)
    for i=1:Nclass
        d_test(j,i) = distance(COVtest(:,:,j),CallClass{i},method_dist);
    end
end

[~,ix] = min(d_test,[],2);
Ytest = labels(ix);
plot(Ytest);
% covariance estimation
disp('done');

%% test 
% figure;
% covWindow=[ones(1,200) zeros(1,200)];
% convulated = conv(d_rest,covWindow);
% 
% convulated = convulated / 200;
% plot(convulated);
% plot(d_rest);
% hold on
% conv(d_rest,covWindow);
% plot(convulated,'r');
% 
% 
% %plot(d_rest>(convulated));
% count = count_sequence(d_rest>convulated(1:length(d_rest)),1);